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Abstract. 

We develop configuration interaction method to treat far from equilibrium quantum 
many-body systems. The concept of embedding (buffer zones between the reservoirs 
and the correlated quantum system) is used to derive an exact master equation for the 
reduced density matrix. To enable the use of field theoretic methods we employ the 
superoperator techniques in Liouville-Fock space. Within the suggested method the 
exact steady-state density matrix is represented as a linear combination all possible 
nonequilibrium quasiparticle excitations built on the appropriate reference state. As an 
example we consider the electron transport through the system with electron-phonon 
interaction. Using approximate (truncated) expansion of the density matrix we obtain 
the linear system of equations for two-quasiparticle amplitudes. Then we compute 
the current and compare the result with other approaches. The current conserving 
property of the method is proved. 
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1. Introduction 

The description of correlated quantum many-body systems far from equilibrium remains 
one of the challenging problems in modern statistical mechanics, quantum field theory, 
nuclear, atomic, molecular and condensed matter physics. [H Ej [3] Most of theoretical 
approaches are based on Keldysh nonequilibrium Green's functions (NEGF) which 
enables perturbative treatment of nonequilibrium many-body systems. [H [3] NEGF also 
allows systematic summation of specific classes of nonequilibrium diagrams, for example, 
random phase approximation |5] or GW theory. [6l [7] Nonperturbative methods, such as, 
numerical renormalization group theory, are also available but often restricted to the 
oversimplified model systems. J8j In this paper we propose a new theoretical method 
- nonequilibrium configuration interaction (NECI), which provides a nonperturbative 
treatment and, in principle, is able to achieve the exact solution of the nonequilibrium 
problem by representing the density matrix as a linear combination of nonequilibrium 
quasiparticle excitations above the vacuum state. 

Let us briefly explain the main ideas of our approach. We start from the 
superoperator representation of the master equation for the density matrix. [9l [101 El 
[12]. A summary of the superoperator formalism relevant to this work is provided 
in Appendix A, while a detailed discussion can be found elsewhere [9l [13]. In 
this representation nonequilibrium dynamics of many-body system is given by the 
Schrodinger-like equation 

^^^\p{t)) = L\p{t)) (1) 

where \p{t)) is the vector in Liouville-Fock space, which corresponds to the density 
matrix, and L is the Liouville superoperator (Liouvillian). In this paper, we focus on 
the nonequilibrium steady-state density matrix, which is a solution of the equation 

L\p) = 0. (2) 

Once the equation is solved, the average of any physical operator A can be computed 
according to 

{A)=TT[Ap] = {I\A\p). (3) 

Here, (/| is the bra- vector in the Liouville-Fock space which corresponds to the 
unit operator /, and the superoperator A is related to an operator A according to 
definition f lA.4p . Thus, our goal is to find the Liouville-Fock space vector which satisfies 



condition 

By the application of Wick theorem in Liouville-Fock space [13] the Liouvillian can 
be brought to the normal ordered form. The diagonalization of the quadratic part in 
terms of nonequilibrium quasiparticles, i.e. in terms of normal modes of the uncorrelated 
part of the Liouvillian, yields: jTH [151 [13] 

L = ^i^nClCn - ^nHi^n) + L' {c\ C,C^ ,c} 

n 

= + L', (4) 
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where c^, and c, c are creation and annihilation superoperators for nonequihbrium 
quasiparticles. The Liouvilhan (jl]) describes a system of nonequihbrium quasiparticles 
with complex spectrum Qn, — ^^e quasiparticle-quasiparticle interaction 

Following the ideas of equilibrium configuration interaction method. [T6] in order 
to describe many-body correlations we need to specify a reference state with respect 
to which the correlations are defined. We choose the reference state as a vacuum 
for nonequihbrium quasiparticles, i.e., c„ Ip*-"-*) = Cn \p^^^) = 0. Then the exact 
nonequihbrium density matrix can be represented a linear combination of the 
vacuum, two nonequihbrium quasiparticle excitations, four nonequihbrium quasiparticle 
excitations and so on. Therefore, NECI expansion for the correlated exact density matrix 
has the following form 

|p) = (l + S)|p(°)), (5) 

where S = J2k contains a sum of the multi-quasiparticle creation superoperators 
which generate /c-quasiparticle excitations above the reference state Ip^^-*), i.e., 5*2 ~ 
~ ^tn^nPp^qy The relative weights (amplitudes) of multi-quasiparticle 

excitations can be found by demanding that the density matrix ([5]) satisfies the steady- 
state condition ([2]), i.e. 

L5|p(°)) = -L'|pW) (6) 

This equation is equivalent to the system of inhomogeneous linear equations for the 
amplitudes of the density matrix expansion. To compute the average of a single-particle 
operator (current, population, etc.) we need to know two-quasiparticle amplitudes. 
However, some truncation are to be made in the density matrix expansion for practical 
purposes. This results in an approximate solution of equation ([6]). 

In this paper, we develop this approach on the specific example of electron transport 
through a quantum system (central region) with electron-phonon interaction. But the 
method can be easily extended to an arbitrary correlated central region. 

The remainder of the paper is organized as follows. In Sec. 2, we introduce the 
Lindblad master equation for an embedded system and its superoperator representation. 
Section 3 presents the nonequihbrium configuration interaction method. We discuss 
two different truncated expansions for the density matrix and apply the theory to 
calculate the steady-state current. In this section, we also proof that NECI method 
is current conserving. The numerical NECI calculations of the steady-state current 
and comparison with other approaches are presented in Sec. 4. Conclusions are 
given in Sec. 5. Appendix A contains a summary of the superoperator formalism. 
Finally, in Appendix B we provide details about NEGF calculations within the first 
Born approximation and the second-order perturbation theory. We use natural units 
throughout the paper: h = ks = \e\ = 1, where — |e| is the electron charge. 
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Figure 1. Schematic illustration of embedding of an open, interacting quantum 
system. Upper part: Each electrode is divided into a macroscopic reservoir and a 
finite buffer zone. The region enclosed by the red dashed line is the embedded system. 
Lower part: The projection of the reservoirs results into the Lindblad master equation 
for the reduced density matrix of the embedded quantum system. Each buffer zone 
energy level Ska is connected to the dissipator and the central region. 



2. Lindblad master equation in Liouville-Fock space 

As shown in Figure [H we consider a paradigmatic model of nonequilibrium many-body 
theory: the interacting quantum system (central region) connected to two noninteracting 
electrodes, left (L) and right (R), maintained with different chemical potentials, /il and 
fiR, and the same temperature T. To be specific, we focus on the electron transport 
problem through one electronic single-particle level with energy Eq coupled linearly to a 
vibrational mode (phonon) of frequency coq (the so called local Holstein model). Thus, 
the central region Hamiltonian is 

Hs = Soa^a + cood^d + fi;a^a((i^ + d), (7) 

where (a) and d^ (d) are electron and phonon creation (annihilation) operators, 
respectively. The left and right electrodes are modeled as semi-infinite tide-binding 
chains of atoms characterized by the hopping matrix element and the on-site energy 
ea (a = L, R). The coupling between the cental region and the end site of a electrode 
is given by the matrix element ta- 

We replace the infinite system by a finite one using open boundary conditions. To 
this end each electrode is partitioned into two parts: the macroscopically large reservoir 
and the finite buffer zone between the central region and the reservoir||| In doing so, 

I The idea of the buffer zone was first proposed by us[14j [151 [13] and then also employed by Ajisaka 
et al.[n] 
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the Hamiltonian of the whole system takes the form 

U = Hs + Hsb + Hb + Hbr + Hr. (8) 
In the energy representation the reservoir and the buffer zone Hamiltonians are diagonal 

Hr = ''*^^eiaa\aO'la-, H B = £ kaol-aO'ka, (9) 



la ka 



where eia denote the continuum single-particle spectra of the left (a = L) and right 
(a = R) reservoir states, while the buffer zones have discrete energy spectrum Ska- The 
buffer-reservoir and central region-buffer couplings have the standard tunneling form: 



Hbr = ^{vikaal^ttka + h.c), (10) 

Ika 

HsB = J2^tkaal^a + h.c). (11) 

ka 

Projecting out the reservoir degrees of freedom from the Liouville-von Neumann 
equation for the total density matrix, we obtain Lindblad master equation for the 
reduced density matrix of the embedded system (central region + finite buffer zones) 

[mils] 

^^ = [H,pit)] + ^Upit). (12) 

Here, H is the Hamiltonian of the embedded system which includes the Lamb shift of 
the buffer zone single-particle levels 

H = Hs + Hsb + Hb + J2 ^kaal^aka, (13) 

ka 

and np(t) is the non-Hermitian dissipator given by the standard Lindblad form 

np(i) = E E (2Wp(^)4., - {Ll,L,^,,p{t)}). (14) 

ka fJ.=l,2 

The operators Lkai and Lka2 are referred to as the Lindblad operators, which represent 
the buffer-reservoir interaction. They have the following form: 



^kal 



with Ffc^i = 7fca(l-/fca), Tka2 = Ikafka- Here fka = [l + e(^*- ^°)/^] ^ and 7^ (Afc„) is 
the imaginary (real) part of the self-energy arising from the buffer-reservoir interaction 

Y.I \vika\'^/{£ka - £la + «0+). 

The Lindblad master equation ( fT2i) describes the time evolution of an open 
embedded quantum system preserving Hermiticity, normalization, and positivity of the 
nonequilibrium density matrix. Open boundary conditions are taken into account by the 
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non- Hermit ian dissipative part, Ilp{t), which represents the influence of the reservoir on 
the buffer zone. The apphed bias potential enters into the master equation via fermionic 
occupation numbers fka which depend on the temperature and the chemical potential 
in the left and right electrodes. We have recently demonstrated that this embedding 
procedure makes the master equation (fT2l) exact in the steady-state regime, if the 
buffer zones are large enough to cure the deficiencies of Born-Markov approximation 
for treating the buffer-reservoir interface. |14[ [T5| [T3] 

We use a superoperator formalism and convert the Lindblad master equation ( |T2l) to 
the non-Hermitian Schroedinger-like equation for the nonequilibrium density matrix. p3| 
[T3] Within the superoperator formalism the density matrix is considered as a vector in 
the Liouville-Fock space and all other operators as superoperators. In Appendix A we 
present a brief summary of the formalism and demonstrate how to convert the Lindblad 
master equation into a superoperator form. The resulting equation is 

|p(t)) = (H-H-tY^ n.«) |p(t)) = L \pit)). (16) 

ka 

Here, the superoperators H = Hs + Hsb + Hb and H = Hs + Hsb + Hb are obtained 
from the Hamiltonian ( IT3l) by replacing ordinary creation and annihilation operators 
a with non-tilde, a^, a, and tilde, a, superoperators, respectively (note that we 
include the Lamb shift into Hb). The dissipators are given by 

Hfca =(Xkal — ^ ka2) {o^ka^ka + (l\a^ka) 

— '^'ii^kalO'kaO'ka + ^ ka2Cl'ka^ka) + '^^ka2- (17) 

The Liouville superoperator L is non-Hermitian because of Ilka- Some other important 
properties of L are discussed in Appendix A. In particular, the relation {I\L = holds, 
where the Liouville-Fock space bra- vector (/| corresponds to the unit operator /. 

We focus our attention on asymptotic (t — )• oo) nonequihbrium steady-state 
situation, where the density matrix \p{t)) does not depend on time. Therefore the 
problem is reduced to the problem of finding the right zero-eigenvalue eigenvector of the 
non-Hermitian finite-dimensional Liouville superoperator 

L\p) = 0. (18) 

Knowing the solution of (fT8l) . we can compute the steady-state current into the central 
region from a buffer zone as 

Ja = {I\Ja\p), (19) 

where 

Ja = -i"^ tka{a\ji - a)aka) (20) 

k 

is the current superoperator. 

In the next section we discuss the development of NECI method to solve 
equation (fT8|) . We also present the application of NECI to compute the steady-state 
current through the central region. 
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3. Configuration interaction method for nonequilibrium density matrix 

In the beginning of the section, we make the important remark on the notation use 
in the rest of the paper: only creation/annihilation super operators written with letters 
a, d (such as for example and are related to each other by the Hermitian 
conjugation; all other creation c\ h\ 7^ and annihilation c, 6, 7 superoperators (as 
well as their tilde-conjugated partners) are "canonically conjugated" to each other, i.e., 
for example, does not mean (c)^ although cc^ ± c^c = 1 (± - bosons/fermions). 



3.1. Nonequilibrium quasiparticles 

Let us introduce nonequilibrium quasiparticle creation and annihilation superoperators, 
which are the basic building blocks for the NECI method. We start by decomposing the 
Liouvillian as 

L - L(°) + L', (21) 
where L^^^ — L^^^ + L^^^ is the LiouviUian for noninteracting electrons and phonons, 

-^ii^ = £o(a^a - a^a) + {Hsb + Hb - t.c.) - i ^ Uka, 

ka 

L^^=uJo{d''d-d)d), (22) 

while L' represents the interaction between them 

L' = K{d^d{d^ + d) - t.c.} (23) 

Hereinafter, the notation 't.c' stands for tilde conjugated superoperators (see the tilde- 
conjugation rules in Appendix A). 

Nonequilibrium quasiparticles are defined by diagonalizing the electronic 
part of L(o): 

L^^=J2(^ndic^-K'^nCn): (24) 
n 

where 1 < n < 2N -\- 1. Using the equation of motion method we find non-unitary (but 
canonical) Bogoliubov transformation which diagonalizes Lg°^ : 

Cn = i^nb + i(Pnb^ + ^{^n,ka^la + i^n,ka^ka)i 
ka 

ci=i^ny^ + ^i^n,kabkai 4 = (Cn)~, 4 = (CnF' (25) 



ka 



where 



P ^a) - ia, b^a, 6^ = {b^)~ b = {b)~ 

ka ~ ^la '^^ka-i ^ka ~ (1 fkaj'^ka '^fka^kai 
^L = (^a)~ hka^ika)- (26) 
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Amplitudes ip, (p are the solution of the following systems of equations 

ka 

and 

^kaVn,ka tkafkai^ 
ka ka 

i^ka ~ ^n)'Pn,ka — ^ka^n = —tkafka'ipn- (28) 



with Eka = £ka + '^ka — ^Ika- By solviug the eigenvalue problem (1271) we also obtain 
the quasiparticle spectrum — fi*. 

Taking into account property ( lA.Tp . we see that (/| is a vacuum state for c]^ and 
c]^ superoperators, 

(/|ct =(/|4 = 0. (29) 

Although superoperators c\^ and c„ (cj^ and c^) are not Hermitian conjugate to 
each other, they obey the fermionic anticommutation relations. Particularly, from 
{c]^, c„} = 1 it follows that ip amphtudes are normalized according to 

(^n)' + $^(^fca)' = 1. (30) 

ka 

Other useful relations between amplitudes ■?/', ({> can derived by making use the 
transformation inverse to ( l25l) (see relations (A4) in our previous paper [ll]) and 
computing the anticommutators. For example 

n 

{^^L} = 0, ^ $^^n,fcaV^n = 0. (31) 



To compute the steady-state curent (]T9l) we need to express the current 
superoperator in terms of quasiparticle superoperators. By doing so and taking into 
account property ( l29i) . we obtain the following expression for the steady-state current 

Ja={I\Ja\p) 

= -2Im ^ tkai'm (y^m,ka + ^ Cfca-^mn) ) (32) 
m,k n 

where F^n = —i{I\cnCm\p) is the two-quasiparticle amplitude. Note, that this is an 
exact expression for the steady-state current. The only problem is to find the unknown 

^ mn- 

In the absence of electron-phonon interaction, i.e., when the nonequilibrium density 
matrix does not contain quasiparticle excitations, the current through the central region 
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is given by the first terms in f l32l) . In what follows we will refer to this current as a 
free-field current Ja'' ■ The correction to the free-field current, AJ^,, is given by the 
second term involving Fmn- In [13], it was shown how to compute Fmn making use 
of the perturbation theory. In the following subsections we demonstrate how to find 
Fmn within the NECI approach using different reference states for the density matrix 
expansion, namely, free-field and coherent reference states. 



3.2. NECI expansion on free- field vacuum as a reference state 

We define the free-field reference state as the steady-state density matrix in the absence 
of electron-phonon interaction, i.e., Ip*-*^^) = |p*-''-')ei |p''°^)ph and 

|pW) = 0. (33) 

The density matrix |p''°^)ei is a vacuum state for nonequilibrium quasiparticles, i.e., 

C„|p(°))el = C„|p(°))el = 0. (34) 

To determine |p*-°-')ph we take into account the possibility that the phonon subsystem 
can contain a certain number of thermally excited vibrational quanta. Let A^^ be the 
number of thermally excited phonons, i.e., 

iV^ = (/|rftJ|p(o)^p^_ (35) 

It is convenient to perform a non-unitary canonical transformation and introduce new 
phonon operators 

^={l + N^)d-Nj\ 

= S -d (36) 

and their tilde conjugated partners 7, j'^ such that |p*-°-')ph is the vacuum state for 7,7 
superoperators, 

7|p^°^)ph = 7|p^°^)ph = 0, (37) 

while (/| is the vacuum for 7^7^ superoperators (see property (1A.7I1 in Appendix). 
Now, the free-field Liouvillian takes the form 

L(°) = 5^(fi„ct - fi:4c„) + oooii^ - 7^7), (38) 

n 

and the vacuum state Ip*-"-*) for annihilation superoperators c„, c^, 7, 7 is the free-field 
reference state obeying condition ( l33l) . The free-field density matrix Ip*-*^^) is normalized 
according to (/|p*^°^) = 1. 

Due to the electron-phonon interaction the exact steady-state density matrix 
contains multi-quasiparticle and multi-phonon excitations above the free-field density 
matrix. We are looking for the exact steady-state density matrix in the form 

|p) = (l + 5)|p(°)), (39) 
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where the operator S is given by the infinite power series of creation superoperators 

oo 

s=Y.a,Ql gl = 4,---4A---4,(7Wr- (40) 

Defined this way, the density matrix is normahzed according to (/ |p) = 1, as (/| is a 
vacuum for quasiparticle and phonon creation superoperators. Since the density matrix 
is tilde-invariant (see the definition in Appendix A), the superoperator S obeys the 
property S = (5')~, i.e. it remains the same if we complex conjugate all and replace 
the non-tilde superoperators by the tilde ones and vice versa. 

We demand that the steady-state density matrix obeys the equation L \p) = 0, 
therefore 

|^L'+ [L(°) +L',5']} = -L'|p(°)). (41) 

This superoperator equation is equivalent to the infinite inhomogeneous system of linear 
equations for amplitudes in the expansion for nonequilibrium density matrix ( l40l) . 
Any truncation in the expansion ( HOi) leads to approximate solution of equation ( l4Ti) . 

Here, we consider the simplest form of S which allows us to calculate the correction 
for the free- field current. Namely, 

S = I Y.{F^n + ^^„7^ + Zl^^^)cl7^^ + ty(7^ + 7^), (42) 

mn 

where Fnm = and is a real number. Inclusion of W and Z terms into the density 
matrix expansion is necessary, since cl^cl"y\ 7^ configurations and their tilde-conjugate 
contribute to the right-hand side of (jUJ) (see the expression for L' below). If we neglect 
them, we observe a homogeneous linear system having a trivial solution. In this sense, 
these terms are correlations inducing terms. 

In order to find equations for the amplitudes F, Z, W we first express L' in terms 
of nonequilibrium quasiparticle superoperators: 

L' =K J2{ K).7^ + L'^nl^ + L^lil + 7)]4cn - t.c.} 

mn 
mn 

-^KY,L^mW - l^YCmCn + /^^^(f - f ). (43) 
mn 

Here the coefficients Lmn are 

Lmn = [(^m " V^m) + A^a;V'm]^n, 

Lmn = [Pm + A^a;^m] ^n, L^^n = ^mi^n, 

L'ml = [i'^m - <fm)fn + N^i'ipmfn " '/'mO] 

Ll^lr = ^mVn - Vm^n. L^^i = ^Pm^n. (44) 
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and n^^^ = {I\d^a\p^^'^) = '^^'ipnfn is the free-field population of the electron level. 
Note, that after normal ordering L' does not involve terms containing annihilation 
superoperators only. Therefore, the condition {I\ L = is fulfilled. 

Substituting S given by (H2l) into (HTj) and demanding equation (14T]) to be fulfilled 
up to terms included into S we derive the system of linear equations for unknown 
amplitudes: 



Fmni^m - ^n) + -^mi [Zin + Z*^] 

i 

i 
i 

Wuo-t^ J2 L'^lFmn = (45) 

mn 

It should be pointed out that the first and the last equations above are exact in the 
sense that the inclusion of other terms into S does not modify these equations. 

3. 3. Coherent reference state for NECI expansion 

Within NECI method we have some flexibility in the choice of the reference state. 
Ideally, we would like to put as much correlations as possible into the reference state 
while maintaining the possibility to define it as a vacuum for some quasiparticle 
annihilation operators. We note that L' given by (l43l) contains the term Kn^^\'^'^ — 7''') 
which is linear in phonon superoperators. It gives us an idea to eliminate this term by 
the non-unitary canonical transformation 

= f , e = 7 + (46) 

and ^"l" = (^^)~, i = The vacuum of the "shifted" phonon operators is the coherent 
state 

|pW) = exp |-«:--(7t +^t) I |p(o)^^ (47) 

normalized according to (/| p*-*-*) = 1. 

Let us consider this coherent state as a reference state for the configuration 
interaction expansion, i.e. 

|p) = (l + 5)|pW), (48) 

where 5" is given by 

00 

^=5^AQ., Q. = <...4,Sj,,...4,(e^r(?)''- (49) 

i=l 
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To distinguish from the NECI expansion on the free-field reference state, we denote 
the present method as NECI*. The advantage of NECI* approach is that it effectively 
includes multi-phonon excitation and de-excitation processes via the exponent (1471) in 
the coherent reference state. 

To fulfill the condition L(l + S) Ip*^*-*) = 0, the superoperator S must obey the 
equation 



where 



and 



|^L'+ [L(°) |pW) = -L'IpW), (50) 

L(°) = ^(n„ct c„ - fi:4c„) + uoie^ - ^^0, (51) 



mn 

+ kJ2{ [L^r^n^^ + L'^ll^ + L^lii + 0]4C„ - t.C.~ 
mn 

— iK ^ [-Z^mn^^ ~ i^nm)*^'^ + ^mlii^ + O] '4nHi 

mn 

-inJ2L^lie-^^)cmCn. (52) 

mn 

If all terms are included into the density matrix expansion, the two methods, NECI 
and NECI*, coincide because both are formally exact. However, when we truncate the 
expansion for the density matrixes, NECI and NECI* give different results. 
Likewise to NECI, we take the operator S in the form 

S = I Y,iFmn + Z^r^S} + ZlJ^^l^^^ + W {i^ + f ). (53) 

mn 

Demanding that equation (150|) is fulfilled up to terms included into S we obtain the 
following system of equations 

2^2^(0) 

+ K ^ [Zin + ^* j] — K ^ {L^ni ) * l^mi + Z*^ 

i i 

— 2kL^^^W = 

Zmni^m " + ^o) y^[-^mi^m ~ (-^niV^mi] 

9^2 (0) 

+ '^Et^J^J^- - ^L^SfFm.] + -^L^^nW = nL^l 
Wu:,-KY,L^lF,r.n = ^. (54) 
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Solving this system we compute the two-quasiparticle amphtudes F^n and, hence, obtain 
the NECI* correction to the free-field current. Here again, the first and the last equations 
are exact, i.e., additions of new terms to S does not modify these equations. 

3.4- NECI is a current- conserving theory 

The truncation of the density matrix expansion introduces approximations into the 
theory. It is important to ensure that observables computed with the truncated density 
matrix still satisfy conservation laws dictated by the symmetries of the underlying 
Hamiltonian. [ini [201 El] In the nonequilibrium case it is particularly important to 
demonstrate that the proposed configuration interaction theory preserves the particle 
number continuity equation in all orders of configurations included into the density 
matrix. To this end, we are going to prove that there is no artificial current leakage 
from the system introduced by the approximations and the current which enters the 
system from the left reservoir is exactly the same as the current which leaves the system 
to the right reservoir: 

Jl + Jr = 0. (55) 

Let us first prove the current conservation in the free-field approximation. Indeed, 
the following equality is true 

^ + = -2Im ^ tka'-PnM'^n 

n,ka 

n,ka 

= Im ^ tkafkai^nMi^n = 0. (56) 
n,ka 

Here, we have used relations (131!) as well as the first equations in (127]) and (!28|l . 

The solution of system ( H5l) (or (1541) ) provides us the two-quasiparticle amplitudes 
Fmn and the NECI (NECI*) correction AJq, to the free-field current (see equation fl32|) ). 
With the help of the first equation in ( 1271) and using the property Fnm = F^n 

AJl + AJr = -2Im ^ tfcaCfca^m-^mn 
mn,ka 

mn 

= -lmJ2^,^rni^m - K)Fmn. (57) 

mn 

Expressing {VLm — VL*^)Fmn from the first exact equation in ( H5l) (or (Ell)) and taking 
into account the explicit expressions for Lmn-, Lmn we find that AJ^^ -|- A J/j = 0. Note 
that this result does not depend on the particular choice of S, i.e. it is valid for any 
truncated NECI (or NECI*) expansion. 
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Thus, the presented configuration interaction theory is current- conserving in all 
orders of nonequilibrium quasiparticle configurations included into the density matrix 
expansion. 

4. Numerical Results 

In our numerical calculations we assume that the applied voltage V symmetrically shifts 
the on-site energies, i.e., e^^R = ±0.5V. Additionally, we assume that the electrodes are 
half filled, i.e., the corresponding left and right chemical potentials are positioned at ei fj. 
We choose the following parameters of the electrodes: Vl = Vji = 2.5, t^ = tji = 1.0 and 
the temperature is T = 0.1. Since the on-site energies e^^R are affected by the apphed 
voltage the left and right electrodes are nonidentical. In the subsequent discussion we 
set the applied voltage V = 1.0. 

For the buffer zones we take a finite chain of N atoms from each electrode. Therefore 
the energy spectrum of each buffer zone is given by 

= e« + 214cos(^^^^J^), k = l,...,N. (58) 

and the tunneling matrix elements in (ITTl) are 

4. = t.y^sin(^). (59) 

The parameter in Lindblad operators (fT5|) is taken to be equal the distance between 
neighbor energy levels in the buffer zones, i.e., jka = ^ka — ^k+ia- In our calculations we 
include = 800 sites into each buffer zone. This size of the buffer zone has been proven 
to give the exact results for the steady state density matrix of the system. [HI [T5] 

Let us say some words concerning numerical solution of the systems of 
equations ( H5|l and (15^ . To be specific we consider the system obtained within NECI 
theory (for the NECI* system we have used the same solution method). The dimension 
of the system is of the order of 4M^, where M = 2A^ + 1. Although the systems is sparse, 
it contains about 12M^ nonzero matrix elements. Assuming that each complex number 
requires 16 bytes of computer memory, we find that for A^ = 800 the total system matrix 
needs about 700 gigabyte for storage. A required memory size can be reduced drastically 
if we contract the indices by introducing the following linear combinations 

fn = Y^ Fmn, = ^(^,nn + Km) (60) 

m m 

and their complex conjugate. Note, that the correction to the free-field current can be 
written as 

AJa = -$^W;fca/n. (61) 

n,k 

Then, it is possible to rewrite the system (145|) as a system for /„, Zn and their complex 
conjugate. The obtained system contains 4M -|- 1 equations and the total number of 




Figure 2. Current through the central region calculated within different approaches 
as a function of the level energy, Eq. 



nonzero elements is about 6M^. To solve the system we have used standard routines 
from Intel MKL Fortran library. 

In figure[2]we compare the second-order perturbation theory (SOFT) and first Born 
approximation (IB A) results for the electron current with the results obtained within 
NECI and NECI* method for different values of the central region parameters. The 
currents are shown as functions of the level energy, Eq, i.e., for different transport regimes 
(resonant /off- resonant). Since the first Born approximation does not guarantee the 
current conservation, we plot the IBA current from the left electrode to the central region 
as well as the IBA current from the central region to the right electrode. The details 
of SOFT and IBA calculations using nonequilibrium Green's functions are presented in 
Appendix B. 

We first consider the case when A^^ = 0, i.e., there is no equilibrium thermally 
excited vibrational quanta. As evident from the figure in this case the first Born 
approximation does not violate the current conservation significantly since the left and 
right currents are close to each other. We see that for a weak coupling {k = 0.5) all 
approaches predict a similar Eq dependency of the current, which reaches a maximum 
value at Eq ~ k^/wq. The peak value of the SOFT and IBA currents slightly exceeds the 
NECI and NECI* ones. When we increase the electron-phonon coupling (k = 1.0) both 
SOFT and NECI currents become unphysical negative in the off-resonant regime when 
the electronic level Eo is below the electrode chemical potentials. In the resonant regime 
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SOPT and NECI currents significantly overestimate the results of other approaches. In 
addition, comparing NECI*-based calculations with other calculations, we can see that 
NECI* approach gives the current peak position at Eq ~ k^/wq. If we now consider 
the case of larger phonon energy {uq = 2.0) we notice that negative SOPT and NECI 
currents in the off-resonant regime disappear and all approaches again demonstrate a 
similar Eq dependency of the current. 

Now we consider results obtained for nonzero number of thermally excited 
equilibrium vibrational quanta. Namely, we assume that = 1.0. In this case the first 
Born approximation clearly reveals its current nonconserving nature. This is illustrated 
in the lower panels of figure [2] where we can see that the left and right IB A currents can 
be very different from each other. Nevertheless, in the weak electron-phonon coupling 
regime {k = 0.5), all approached predict qualitatively similar dependence of the current 
upon the electronic level Eq. Comparing the result for k = 0.5, N^^ = 1.0 with those for 
A = 0.5, iV^ = we notice that inclusion of thermally excited equilibrium vibrational 
quanta into consideration produces slightly broadened current peak with a reduced 
amplitude. When the coupling constant increases {n = 1.0) SOPT and NECI approaches 
again show unphysical negative current in the off-resonant regime when Eq is below the 
electrode chemical potentials. Moreover, for n = 1.0 the electron-phonon coupling splits 
the NECI* current peak, while NECI approach gives only one peak. For a larger phonon 
energy, uo = 2.0, NECI* predicts three pronounced peaks, while NECI gives only two 
peaks. In addition, SOPT and IBA show unphysically negative current for certain 
values of Eq. 

From the above consideretion it is evident that the configuration interaction method 
built on the coherent reference state is preferable when the electron-phonon coupling is 
large or when there are thermally excited vibrational quanta in the system. This result 
is not surprising since, contrary to NECI, NECI* method accounts for multi-photon 
excitations and de-excitation processes which are important in these regimes. 

5. Conclusions 

We developed nonequilibrium configuration interaction method, which formally gives 
the exact solution of out-of-equilibrium correlated many-body problem. Our approach is 
based on a superoperator representation of the Lindblad master equation for the reduced 
density matrix of the embedded quantum system. It was shown that the steady-state 
density matrix can be decomposed in Liouville-Fock space in terms of nonequilibrium 
multi-quasiparticle excitations above the reference vacuum state. The amplitudes of 
these excitations provide a measure of the many-body nonequilibrium correlations. 

The theory was applied to study the inelastic electron transport through the system 
with electron-phonon interaction. To compute the current we used truncated expansion 
of the steady-state density matrix. Two different reference states were considered: free- 
field vacuum and coherent state. It was proved that both approximations are current 
conserving in all orders of the density matrix expansion. The current through the 
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system was computed for different model parameters and compared with the second- 
order perturbation theory and the first Born approximation results. It was shown that 
the configuration interaction method based on the coherent reference state is superior 
to the other approaches when the electron-phonon coupling is large or when there are 
thermally excited vibrational quanta in the system. 

The method can be readily extended to dynamical nonequilibrium case by making 
the amplitudes in the density matrix expansion time-dependent functions. 

Appendix A. Brief overview of the superoperator formalism 

Here we briefly review the formalism of superoperators acting in the Liouville-Fock 
space, while the detailed discussion is given in [9l [13]. We suppose that the Fock space 
of the system under consideration is spanned by the complete orthonormal set of basis 
vectors |m) = |mi,m2, . . .) which are eigenvectors of the particle number operator: 

a\ak\m) = mk\m) , ^ |m)(m| = /, {n\m) = 5nm- (A.l) 

m 

The operators in the Fock space form themselves a linear vector space called the 
Liouville-Fock space. The set of vectors \mn) = \m){n\ constitutes a orthonormal basis 
in the Liouville-Fock space. Thus, every Fock space operator A = J2mn \m){n\ can 
be considered as a Liouville-Fock space ket-vector \A) = Ylimn^mn \mn). The adjoint 
operator is represented by the bra- vector {A\ . 

The scalar product in the Liouville-Fock space is defined as (A1IA2) = 1:i{A\A2). 
In particular, the scalar product of a vector \A) with (/| is equivalent to the trace 
operation in the Fock space, {I\A) = Tt{A). 

Now we introduce creation and annihilation superoperators. Fermion creation and 
annihilation superoperators are defined as 

dk \mn) = at \m)(n\ , Ht \mn) = «(— 1)^ \m){n\ a|. 
d\ \mn) = a\ \m){n\ , a\ \mn) = z(— 1)^ \m){n\ a^, (A.2) 

where /i = X]fc(^fc + ^fe)- For bosonic creation and annihilation superoperators we drop 
the phase factor z(— 1)^. So defined superoperators satisfy the same (anti) commutation 
relations as their Fock space counterparts. Additionally, the basis vectors of the 
Liouville-Fock space are super-particle number eigenvectors, i.e., 

dldk \mn) = \mn), a^a^ \mn) = ut \ran). (A. 3) 

For an operator A = A{a\ a) we formally define two superoperators 

A = A{a\d), A = A*{Z\a) (A.4) 

and refer to them as non-tilde and tilde superoperators, respectively. The connection 
between non-tilde and tilde superoperators is given by the "tilde conjugation rules" 

(ciii + C2A2y= clAi + CIA2, 

{A,A2r=A,A2, {Ay=A, {A^r={A)K (A.5) 
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With the help of superoperators an arbitrary Liouville-Fock space vector can be 
represented as 

\A) = A\I) = aAA^\I). (A.6) 

Hereinafter, the phase a a = —i if A is a fermionic operator and (Ja = +1 if ^ is a 
bosonic operator. We also define a Liouville-Fock space vector tilde conjugated to a 
given one, 1^)^^= A\I) = o\A^ Therefore, if A is a Hermitian bosonic operator then 
\A) is tilde-invariant, i.e., \A) = \Ay. The examples of tilde-invariant vectors are |J) 
and the density matrix \p) = p\I). 

Considering the adjoint of flA.6l) and assuming that A = a or we find 

(/| (a^ - <a) = (/| (at - a„a) = 0. (A.7) 

This gives us the idea introduce superoperator = a) — a* a and its tilde conjugate, 
b'^ = a^ — aaa, which annihilate the bra-vector (/| . 

For the product of operators the following relation is valid 

\A,A2) = A^\A2) = TAl\A,), (A.8) 

where r = z if both Ai and A2 are fermionic and r = cr^j otherwise. With the help of 
relations flA.SP the average of an operator A in the state with the density matrix p can 



be calculated as the matrix element of the corresponding superoperator A sandwiched 
between (J| and \p): 

{A) = Tr(Ap) = {I\Ap) = {I\A\p). (A.9) 

Let us now consider the Lindblad master equation f|T2|) . In the Liouville-Fock space 
this equation takes the form 

z| \p{t)) = \Hp{t)) - \p{t)H) + z \Up{t)). (A.IO) 

Applying raltions flA.81) . i.e., taking into account that \Hp(t)) = H\p(t)), \p{t)H) = 
H \p{t)), \akap{t)alj = aka \p(t)alj = iakaO-ka |p(i)), etc., we can rewrite this equation 
in the Schrodineger-like form f[T^ . The time evolution of the density matrix is governed 
by the non-Hermitian Liouville superoperator L. Two important properties of L should 
be noted: 1. Since the density matrix is tilde invariant, then {LJ = —L; 2. Taking 
the time derivative of (/|p(t)) = 1 we find (/| L = 0, i.e., (/| is the left zero-eigenvalue 
eigenvector of L. 



Appendix B. Details of NEGF calculations 

In this Appendix we briefly describe the details of transport calculations based on 
nonequilibrium Green's functions relevant to our discussion. We first define four kind of 
interacting Green's functions: G^'°'{u)) and G^{u)). The retarded and advanced Green's 
functions are given as the solutions of Dyson's equation 

= Gl + G^^^G', G" = (G"")*, (B.l) 
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where G'q = [u — Eq — T,^^ — is the noninteracting retarded Green's function. The 

lesser and greater Green's functions satisfy the Keldysh equation 

= G^(Sf + S| + J:fjG\ (B.2) 

The interaction self-energies S?^^ account for the electron-phonon coupling in the central 
region, while the electrode self-energies, describe the effects of the coupling between 
the central region and the electrodes 

S<(u;) = +2zImK(u;)](/,H-l), 

S>(u;) = +2zImK(u;)]/„H. (B.3) 

Here, fa is the Fermi-Dirac distribution function with chemical potential fia and 
temperature T. 

Once the full interacting Green's functions are defined, we can write explicit 
expressions for physical quantities of interest. For the net electric current passing 
through the a electrode to the central region we adopt the Meir-Wingreeen transport 
formula ^22j 

Ja = J ^{S<(a;)G>(a;) - S>(a;)G<(a;)}, (B.4) 
and for the electron level nonequilibrium population we use 

nei= / p.G<ico). (B.5) 

Neglecting the electron-phonon interaction, i.e., assuming = 0, we obtain the free- 
field current, Ja \ and electron populations, n^^\ 

The exact electron-phonon interaction self-energies in Green's functions (1B.2P 
and (IB.ip contain all possible diagrams that satisfy Feynman rules. In this work, we 
calculate S?^*^ from the lowest order self-energy diagrams referred to as the Hartree 
and Fock diagrams, T,^[ = T,^^ -\- S^''^. This approach is known as the first Born 
approximation (IBA). Using the standard rules of Feynman's diagrams and the Langreth 
theorem for analytical continuation we find the Hartree and Fock self-energies 



(c) = - ^K'Dtico = 0) 1 ^ G<(a;') 



+ D^,{u-u')[G<{u') + Gl{u')]]. (B.6) 
With the usual definition for the bare phonon Green's functions 

Df{uj) = -2m[NJ{uj T ujq) + (A^^ + l)6{u ± uq)] 

DUuj) = -, (B.7) 
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where N^^ is the equihbrium occupation number of the phonon mode with frequency ooq. 
We see that = because Dq{uj = 0) = unless Uq = 0. Since Dq{uj = 0) = —2/ujq, 
we have = —2K'^n^^^ /coq. 

Substituting self-energies (1B.6P into Dyson equation (]B.1|) and Keldysh 



equation ( 1B.2I) we find the lesser and greater Green's functions within the first Born 
approximation. Although the IBA does not generally guarantee current conservation, 
computationally it is not as demanding as the self-consistent Born approximation when 
the full interacting Green's functions are used in the definition of the self-energies ( IB. 61) . 

In the case of weak electron-phonon coupling the second-order perturbation theory 
(SOFT) may be a good approximation. In this approximation the lesser and greater 
Green's functions obtained within IBA are expanded with respect to electron-phonon 
coupling up to terms 

+ 2Gf(a;)Re[GSMS[„,(u;)]. (B.8) 

Substituting this expression into ( IB. 41) we evaluate the current within the second-order 
perturbation theory. It should be noted that this approximation satisfies the current 
conservation condition. 

In the end of this Appendix we present explicit expression for the electrode retarded 
self-energy. For the semi-infinite one-dimensional chain with on-site energy and 
hopping parameter Va coupled to the central region with hoping matrix element ta, 
the retarded self-energy ^^(co') = Aa{uj) — 0.5iTa{uj) can be written in an analytical 
form as 



AJuj) 




X, \x\ < 1 

X — sgn{x)y/x'^ — 1, > 1, 



T^{u) = Fo,e(l - \x\)Vl^, (B.9) 
where x = {u - e„)/(2|K.|), Foa = 2tl/\Va\. 
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